#ifndef InteractionCalculator_Lin_H
#define InteractionCalculator_Lin_H

#include <vector>
#include <set>
#include <map>
#include <cassert>
#include <iostream>
#include <fstream>
#include <sstream>
#include <utility>
#include <ctime>

#include <math.h>
#include <stdio.h>

#include "LBWrapper.hpp"
#include "Sphere.h"
#include "../sMathLib/Vector3/Vector3.hpp"
#include "Box3D.h"
#include "Properties.h"

#include "volumeFractionField.h"

#include "ImmersedObjectWriter.h"

#include "ImmersedObject.hpp"

#include "../../../../vedo/src/NJR/Interfaces/vector3d.h"
#include "MpiFacade.h"

#include "olb3D.h"
#ifndef OLB_PRECOMPILED // Unless precompiled version is used,
  #include "olb3D.hh"   // include full template code
#endif

using namespace std;
using namespace olb;
//using namespace ibm;

//typedef double DecType;
//typedef double T;

typedef unsigned Index;

template<typename T, template<typename U> class Lattice>
struct ClearPhiF : public WriteCellFunctional<T,Lattice>
{
  virtual void apply(Cell<T,Lattice>& cell) const
	{	
		//cell.clearPhi();
	}
};

template<typename DecType>
class InteractionCalculator_Lin
{
public:
	InteractionCalculator_Lin();
	
	InteractionCalculator_Lin(DecType tau);

	~InteractionCalculator_Lin();


	template<typename T, template<typename U> class Lattice>
	void calculateInteraction(BlockStructure3D<T, Lattice>& lattice, T Nu, int iT);

	template<typename T, template<typename U> class Lattice>
	void parallelCalculateInteraction(MultiBlockLattice3D<T, Lattice>& lattice, T Nu, int iT);
								
	void clearAll();
	
	void addDO(NJRvector3d particleCenter, NJRvector3d particleVelocity, NJRvector3d particleAV, DecType r);
	
	unsigned getNumberOfDO() const;
	
	const NJRvector3d& getForces(Index i) const;
	
	const NJRvector3d& getTorque(Index i) const;
	
	void assignBoundaries(int xLhs, int xRhs, int yLhs, int yRhs, int zLhs, int zRhs);
	
	template<typename T, template<typename U> class Lattice>
	void assignComputationalDomain( MultiBlockLattice3D<T, Lattice> &lattice );
	
	void setTau(DecType tau)
	{
		this->tau_ = tau;
	}

	sMathLib::Vector3& getNormalStress(unsigned sphereIndex) const;

	sMathLib::Vector3& getShearStress(unsigned sphereIndex) const;

	void recordFluidEffect(int iT, int recordStep, std::vector<DecType>& scalingFactor);

	template<typename T, template<typename U> class Lattice>
	void createVolumeFractionField(MultiBlockLattice3D<T, Lattice>& lattice);

private:
	
	//for DO
	std::vector<NJRvector3d> linearForce;
	std::vector<NJRvector3d> torque;

	std::vector<DecType> sphere_forces;
	std::vector<DecType> sphere_torque;

	std::vector<DecType> normal_stress;
	std::vector<DecType> shear_stress;

	std::vector<DecType> sum_forces;
	std::vector<DecType> sum_torque;

	std::vector<DecType> sum_normal;
	std::vector<DecType> sum_shear;

		
	std::vector< ibm::Sphere<DecType, ibm::ImmersedObject> > _sphereString;
	
	int LX, LY, LZ;
	
	int counter;
	
	DecType tau_;

	sMathLib::VectorX<int, 3> _computationalDomain;
	Box3D<DecType> _domainBox;

private:

	std::ofstream stressProfile;

	ibm::VolumeFractionField< std::map<int, DecType> >* volumeFraction;

private: /// Parallel version
	
	template<typename T, template<typename U> class Lattice>
	void calculateVfAndModifyPopulation(MultiBlockLattice3D<T, Lattice>& lattice, int iT);
	
	template<typename T, template<typename U> class Lattice>
	void calculateForce(MultiBlockLattice3D<T, Lattice> &lattice, T Nu, T drift_in, T drift_out, int iT);

private: /// Seirial version
	
	template<typename T, template<typename U> class Lattice>
	void calculatePhiAndModifyPopulation(BlockStructure3D<T, Lattice> &lattice, double drift, int iT);
	
	template<typename T, template<typename U> class Lattice>
	void calculateForce(BlockStructure3D<T, Lattice> &lattice, double Nu, double drift, int iT);

};



#endif /*InteractionCalculator_Lin_H*/
